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Unsuccessful cosmology with Modified Gravity Models 

Antonio De Felice*, Mark Hindmarsh^ 
Department of Physics & Astronomy, University of Sussex, 
Brighton BNl 9QH, United Kingdom. 



A class of Modified Gravity Models, consisting of inverse powers of linear combi- 
nation of quadratic curvature invariants, is studied in the full parameter space. We 



o 

. find that singularity-free cosmological solutions, interpolating between an almost- 



Friedmann universe at Big Bang Nucleosynthesis and an accelerating universe today, 
exist only in a restricted parameter space. Furthermore, for all parameters of the 
models, there is an unstable scalar mode of the gravitational field. Therefore we 



\ conclude that this class of Modified Gravity Models is not viable. 
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I. INTRODUCTION 



^ \ One of the most puzzling problems in contemporary physics is the the accelerating expan- 

Xf^ • sion of the universe [1-10]. Standard General Relativity (GR) can accommodate acceleration 

: 

. through a cosmological constant but gives no explanation of its size. 

rn 

in 
o 
i> 
o 



The problem of such a universal constant is difficult to attack with only GR and the 
Standard Model of Particle Physics (SM). This leads to two possibilities, of which one is to 
try to see if extensions to the Standard Model motivated by other physics can lead, more 
or less naturally, to the solution of this problem. Failing that, one can introduce a pure 
phenomenological model, more or less inspired by some symmetry principles, to try to fit 
the acceleration data with the fewest free parameters. 

Among these models, the most developed one is quintessence [11-25], that is a mini- 
mally coupled scalar field with a properly chosen potential. However a completely different 
approach has recently emerged, shifting the acceleration behaviour from the matter sector 
to the gravitational one. In other words, the acceleration, according to these models, is a 
purely gravitational effect, induced by changing the Einstein-Hilbert action for gravity. 
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These modifications must become important at late times in order to explain the fact 
that the universe has started to accelerate only recently, and so much attention has been 
focused on actions with inverse powers of curvature invariants, mainly squares of the Ricci 
scalar i?, the Ricci tensor i?^,^, and the Ricmann tensor R^^ppa [26-37]. These include theories 
where the Einstein-Hilbert action is supplemented by a function of the Ricci scalar alone, 
/(i?), and also /(i?, P, Q) where P and Q are the squares of the Ricci and Riemann tensors 
respectively. 

The first class of theories, in general, can be mapped to a scalar tensor theory [38, 39], 
and are subject to constraints from both solar system tests of gravity ([40-42, 42, 44]) and 
cosmology ([45-49]). The second kind of theories seem to be inspired by the kind of higher 
order terms induced by radiative corrections [50]. It has been noted that on Minkowski 
and de Sitter background, these theories would have ghosts unless P and Q appear in the 
Gauss-Bonnet (GB) combination, that is f — f{R,Q — 4P) [51-60]. While allowing the 
possibility of an accelerating universe today, albeit at the cost of an unexplained small 
parameter, these theories are afflicted by ghosts, instabilities and superluminal modes in the 
accelerating background for a large part of the parameter space [61]. 

If the late-time behaviour constraints future instabilities, in this paper we address further 
bounds based on cosmological behaviour in the past history of the universe. We find that 
the story is no more encouraging. 

In half the parameter space, the dynamical equations, expressed in terms of linear com- 
binations of and Q — AP, can be shown to possess a separatrix corresponding to a scale 
factor a{t) oc t^', with 1/2 < p < 1, which cannot be crossed. This precludes good cos- 
mological behaviour for these models, since starting from a radiation-dominated universe at 
nucleosynthesis, with p ^ 1/2, it becomes impossible for the universe to reach an accelerating 
phase. 

For the remaining part of the parameter space, we have performed a detailed numerical 
analysis in order to find a background which could mimic GR from Big-Bang Nucleosynthe- 
sis (BBN) up to present. Even though it is impractical to integrate the equations of motion 
for such a long interval of time as the modified Friedmann equation is extremely stiff, we 
establish that GR-like initial conditions can be imposed and integrated for small time in- 
tervals. However, we also find that models with inverse powers of quadratic invariants, for 
any value of the power-exponent, possess in general classical instabilities during radiation 



3 



domination. Therefore, we conclude that this class of models cannot be a viable explanation 
for late-time accelerated expansion. 



II. THE MODELS 



We will consider models defined by the following action (see [27]) 

/r B ,,4n+2 

/O ,,4n+2 1 



where b = ai/ag - 1, (ag 0), = l/(87rG'), = sign(//), i^^g = i?^ - 4P + g, and 



(1) 
(2) 
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Q — Raf3jS R"^"'^ , P — Rafj R"' 

For our purposes, it is more convenient to introduce auxiliary scalar fields A and 
rewrite the action as 



(3) 
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0n+l 0n + 0n+l GB^ 



and 



(4) 



where C = 0^ ii'^'"'^'^ / d''^. We will also find useful the form 



S^Ml j d'x [(i + x(A, 0)) i? - U{\, 0) + RIb] , 



(5) 



where 



u = 

e = 



2bnCX 

c 



0" 

0n+l 



on — + n + 1 



6 + 



n + 1 
n P 



(6) 
(7) 
(8) 



With this choice it is clear that the equations of motion, for any background, become of 
second order for the two scalar fields A, 4> and for g^i,, namely 



A = i?. 



bR^ + i?GB ) 



(9) 
(10) 
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and 

(l + X) Ra0 - VaVpX + - 2 V c^V + 2 Q^p R 

-4 V, -i9ap[{l + x) R-U]=A7rT^p . (11) 

Eqs. (9) and (10) are second order differential equations for the metric tensor, where the 
last one involves both second order derivatives in all the dynamical fields, A, 0,(7^1,. We 
will suppose that TaiB is represented by a perfect fluid with two components, radiation and 
coUisionless matter, and that the spacetimc interval has the Friedmann-Lemaitre- Robertson- 
Walker (FLRW) form ds'^ = -dt'^ + a?{t){dx^ + dy'^ + dz^). 



III. THE EQUATIONS OF MOTION ON FLRW 

It is also convenient to perform a change of the time variable, so that 

dt^^, (12) 

and T is dimensionless. Denoting dimensionless quantities whose prime derivatives use r 
with a bar, we have H ^ HqH, ( ^ and 

R — Hq R , X — Hq a , (13) 

RIb-H^RIb: <t> = H',4>. (14) 

Having rescaled variables so that they are dimensionless, we may remove the bars on the 
understanding that H = a'^da/dr, R = Q{H'^ -\-a'^d^a/ dr"^) and Rq^ = 2AH'^{a~^d^a/ dr'^). 
If we parameterize the FLRW metric in the following way [37] 

^ dN^ + al e^^ {dx^ + dy'^ + dz") , (15) 

where H — (3 and N — ln(a/ao), so that w(0) = — ln/3 (defining the dimensionless 
parameter /3), then one has 

R = 6/3'e'"K + 2] (16) 

Rl^ = 24/?^e^"K + l] , (17) 
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where u' — du/dN. Prom the equations for A (or x) and one can find a relation among 
the dynamical variables, that is 

= 6A' + 4/32e'"(A-6/32e2") . (18) 

This means that we need only two independent dynamical variables, which we can freely 
choose the most convenient ones. Now, we have that 

= :S ' (19) 

then u' < 0, assuming that the universe docs not superacceleratc, i.e. assuming H < 0. It is 
interesting to note that Eq. (18) can be solved for /3^e^", obtaining 



^.,^,,,„^A±VA!+|ME1. (20) 

The quantity inside the square root is always positive, as can be easily verified: 

A" + 6(6A2 - 0) = i?" - 6i?^B = 36 e^" u'^>0 . (21) 

Furthermore, in the past, the universe was decelerating and so Rq^ < 0. Therefore we 
take the positive root in Eq. (20), from which we learn the dependence of H on the other 
variables A and 0. Let us consider Eq. (16) for A, which can be written as 

A = i? = 3/32— - + A + VA2 + 6(6A2-0) . (22) 
aiV 

By differentiating Eq. (20) we find 

ly + i (1 + 6^)^^^ -I +Va^ + 6(6P-0)=O. (23) 

' ' VA^ + 6(6A2 - 0) ^ + 6(6A2 - 0) v ^ v^; \ j 

Prom Eq. (11) we can derive the modified Priedmann equation , which is 

8/3V"^' + 2/32e2"(x' + x)-|t^=-^-/3'e'" , (24) 

Pcr,0 

where ^' = (f)', and x' = 26 A (f>' + 2b^ A'. 

Eqs. (23) and (24) can be combined to give a system of two first-order differential equa- 
tions from in the variable A and 0: 

AnX' + Ai2(l>' = Bi (25) 
A21 A' + A22 0' = , (26) 
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where 



A^^^l + ^MB^ (27) 

Ai2 = ^ = (28) 

Si = -4yA2 + 6(6A2 - 0) (29) 

A21 = I fee [A + v/A2 + 6(6A2-0)] (30) 

^22 = II [A + VA^ + 6 (6 A2 - 0)]' + i 6 [a + VA2 + 6 (6 A2 - 0)] (31) 



+ i t/ - I [A + VA2 + 6(6A2-0)] [1 + 46^ A] . (32) 

Therefore one has 

, , ^99 -Bl — ^19 -S9 . . 

A' = ^ (33) 

,/ All B2 - A21 Bi 

^ = > (34) 

where A = An A22 — A\2 A2\. These arc the dynamical equations to solve numerically. The 
choice of these dynamical variables seem to be a natural one, since the two variables have a 
clear geometrical and covariant meaning. 



IV. ANALYSIS OF PARAMETER SPACE 

The action has three free parameters C, ^, 'n- Prom the definition of the variable one 
can see that 

= 12/3^ e^" [36 m'' + 2(1 + 66) («' + 1)] , (35) 
and li' = — (1 + g), where q is the deceleration parameter, so 

= 12//^ [36?^- 2 (1 + 36)^ + 36] . (36) 

Therefore may vanish for particular values of q. It should be noted that if 6 < —1/6, 
cannot vanish, and remains negative for all values of q and therefore at all times. The 
dependence on 6 allows us to divide the parameter space into three regions. 
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A. b>0 



The variable vanishes when 

?i,2 = 3^ . (37) 

provided that 6 > — |, and 6 7^ 0. It should be noted that for the particular case 6 = 0, 
there is only one real solution, q = 0, and (p = Rq-q. For b > —1/6 there arc two distinct 
real solutions. When 4> approaches zero, the equation for 4> can be approximated by 

+ 0(0^) = ,+0(0=^) (38) 

where we have used the fact that A does not vanish at the points where does. The solution 
near = is ^ = Ae"^^ ^ from which it is clear that never crosses zero for finite N . Thus 
the line = is a separatrix in the (A, (j)) plane. 

Suppose we start at = 0. Then a flow is deflned by the equation for A, which must 
have constant q = q^. A constant deceleration parameter implies 

a = d. . (39) 

This has two possible solutions 

a = a, e^* if = -1 , (40) 

a = a*[l + //.(l + g.)(i-t.)](Vi+9«) ifg, ^-1, (41) 

that is either an exponential or a power-law a oc with an exponent defined as p = p* = 
(1 + g*)~^. Thus we learn that the power law expansion parameter p can never cross the 
value p*. Hence for a universe which moves from a radiation-dominated phase p — 1/2 to 
an accelerating phase p > 1 we must not have p* in the interval (|, 1]. 

Therefore we can rule out all models of the type Eq. (2) with 6 > 0, since 1/2 < 
1/(1 + qi) < 1, i.e. there is a separatrix between radiation domination and the time of 
acceleration domination. 



B. 6<-i 



For 6 < I we cannot use the argument of the previous section as the = separatrix 
corresponds to power law expansions outside the range | < p < 1. 
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Let us instead define a new variable k — On dimensional arguments, we might 

expect that k ~ On the other hand let us consider what happens to the equations of 

motion introducing this ansatz. Since 

R^6H'{l-q), (42) 

and q < 1 during and after radiation domination, then 9x — 1 and we have 

^„ = l + ii^ (43) 

A22 = ^ (1 + [1 + 6^ + 7«] (44) 

>li2 = -T^ (45) 
A 7k 

A21 = i6eA(l + 7K) , (46) 

where 7^ = Vl + 66 — 6k,. The definition of k implies that k, < 1/6 + 6. Therefore one has 
for the variable A defined in Eqs. (33, 34) 

^ = 4^^[l86«:-(n+l)(7« + l + 66)'] . (47) 

It is important to check where the zeros of A are placed, as the dynamical equations (33,34) 
become singular at A = 0. This means that the derivatives A' and 0' diverge, and hence a 
divergent derivative of R and Rq^- 

Denoting the value of k, at which A vanishes as k,^, defining ip = (p — k.^ X^, and expanding 
A about (giving A (^kA)^^ V'/AJ + . . . ) one has 

ijj 

assuming that neither F nor (^kA)^^ vanish, which has solutions i/)'^ oc N — N^. This says 
that the solution cannot cross the singularity, because either N > or N < N^. We 
relegate the exceptional cases in which our assumptions are not valid to Appendix A. 
The equation A = has two possible solutions for k, which solve the quadratic equation 

3^-^ + 1^ K^ + 2{l + 6b) 1 + 66- (1 + 36) ^^^^ + 1^ k + 36^(1 + 66)^ = . (49) 

If 6 = — I (n + 1) then k — — — §6^, otherwise there are, in general, two solutions for 
A = 0, namely 

^. = 'j^^^^[n-s^ + V^ITTW^TT)\ (50) 

--= '(lt36+t)?' t"-^^-^ -'^^^^)P"^^] - 
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These solutions are real and negative if6<— |. If6 = — |, then kx — K2 — 0. It should be 
noted that 

Ki<K2<\ + h . (52) 

It is useful to determine the relation between q and n. Using the definitions of A and (j) 
we find 

_ 36g^-2(36+l)g + 36 

%^Tp ■ ^^^^ 

For 6 < — |, K is always negative. According to the evolution of our universe from BBN 
to the point of change of sign for the acceleration (some time in our past), the variable q 
smoothly evolves from 1 to 0, that is 1 > g > 0. In this interval k monotonically decreases 
as q increases, and smoothly changes in the interval — oo < k < b. Therefore if ki < b, 
the universe should have evolved from a radiation dominated universe at BBN to a state 
with a singularity in a derivative of a curvature invariant. Of course this scenario is not 
consistent with our previous analysis as the universe would not be able to cross in general 
such singularity. This happens when 



Ki < b or K.2 <b . 



(54) 



The inequality ki <b can be recast in the following form 

-{n +!)(! + 66)V-(l + 66)(2n+l) > -(1 + 2n)[9b'^ - (n + l){3b + 1)] , (55) 

therefore if 96^ — (n + l)(36+l) > 0, then this is always verified and there is a past-singularity 
for all values of b such that 

n + 1- y/{n + l){n + 5) 



b < 



6 



In the case of 



96^ - (n + 1)(36 + 1) < , or 



n + 1 - ^y{n + l){n + 5) 



<b< — 
6 - - 6 

then taking the square we find that the Ki < 6 is equivalent to saying 

2(n+ 1) - ^2{n+l) 



(56) 



(57) 



{l + 3b + ny 



b + 



2(n+l) + v/2(n+l) 



3(2n + 1) 
This is verified by the interval 

n+1- ^y{n+l){n + 5) 
6 



b + 



3(2n+ 1) 



< . 



(58) 



< b < 



-2(n+ 1) + V2(n+ 1) 
3(2n + l) 



(59) 
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Therefore, in total we have that there is a pathological behaviour for 

2(n + l) + V2(^ (60) 
3(2n + l) ^ ' 

It should be noted that — | < & < — |. The inequality K2 < b gives a weaker bound on b, 
with the past-singularity interval b < — [2(n + 1) + •\/2(n + l)]/[3(2n + 1)] < b. 



C. 6<6<0 

In this region, the differential equation possess no singular points, and therefore one 
can try to solve the equations from BBN up to present, to check that there is a viable 
cosmological solution. We will assume that above a certain redshift {N ^ —20, for BBN) 
the dynamics of the universe follows that predicted by GR very closely, so that in Eq. (24) 
we have 

_/32g2«^Q (61) 



Pcr,0 

In this way wc know how to impose initial conditions for these cosmological models, and we 
will require that they be free of ghosts and have a positive square of the propagation speed 
at the initial "GR-like" time. 

As H'^ — e^" and p oc exp(— 4A^), the scale factor a must expand as a power of r which 
is very close to 1/2: writing a oc t5+^ (recalling that t — Hot) we find 

A = 6e/r^ = -3/2T^ (62) 

and hence that k = 0/A^ <^ 0. Indeed, we checked that for values of N close to BBN, one 
can choose initial conditions for which we have a GR-like evolution for which —8 
and A'/A —4. This can be done as follows 



A, = 6 



<f>,^bX^ + 2AHf^, (64) 



(63) 



where 



Hf = n^o (1 + Zif + ^0 (1 + Zi)^ + 5 (65) 

S-lnmo{l + Zif-nro{l + Zif. (66) 



a. 



11 



In the previous 2 equations, the Priedmann equation and the second Einstein equation for 
the initial conditions have been perturbed by a small parameter 5, equivalent to a tiny 
cosmological constant. Therefore, at a given redshift, the quantities and X'JXi become 
functions of 6. Then we numerically found that particular 6 which gives X'i/Xi = —4 with 
a chosen accuracy. Automatically wc also found that ^y^j = —8, satisfying the property 
that K — is approximately constant ~ —10^°). In order to find such a S, and more 
in general to solve numerically the dynamical equations, one needs multiprecision analysis, 
as must cancel the usual matter (radiation plus dark matter) up to terms of order of 
U <^ H'^. Along the same lines, the full solution of the differential equations requires time- 
steps so tiny that the cvohition (from BBN or radiation-domination to recombination time) 
of these backgrounds becomes impractical (at least at these redshifts). Even though we 
were not able to achieve a full time-evolution for the solution, however we could set up GR 
initial-conditions for these models and check that background was still GR-hke for at least 
a few thousand time-steps. 

On the other hand there are further conditions that should be fulfilled by the propagator 
of the gravitational field: there should be no ghosts and the square of the propagation speed 
should be positive [61]. These conditions can be written as follows 

1. Absence of ghost-behaviour for the spin-2 graviton 

l + 46^A + 8^'>0 (67) 

2. Absence of classical instabilities 

Co > where ct, — r < 1. (68) 

l + 4:bCX + 8HC 

3. Absence of classical instabilities for the spin-0 modes of the metric 

32 8 

cl>Q where d = 1 + - ^^{^ - i H), (69) 

where 

Qi^Ab(iX + ^X) + 8iH^ and Q2^l + ^b^X + 8H^ . (70) 

We now have sufficient information to study the propagator constraints (see [61] for 
details) as functions of the initial conditions imposed at redshift Zi. 
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Noting that A ~ HX, ^ ~ —H^, and that A ^ if^, we find that 

Qi^8^H\ (71) 

whereas 

^2 ~ 1 , (72) 

as ^A and H'^^ both tend to zero as r — > 0. Therefore, since ^ oc H"^^, it is easy to see that 

2 , , 2 . 32 ^ii" ^ 8 5 
Co 1 and Cn fa 1 + A = 1 - - = — . (73) 

The presence of this classical instability makes the evolution of these models inconsistent 
with the presence of an early radiation-dominated phase at the era of BBN, and there- 
fore they cannot be accepted as viable models, even if the background could evolve from 
radiation-domination to matter-domination in a way consistent with today's data. 

The result is that the parameter space of the initial conditions at redshift Zi, which satisfy 
all the no-ghosts/instabilities constraints is empty. In particular the scalar modes have an 
imaginary speed of propagation. This way avoids the formidable problem of solving this 
system of stiff differential equations. 



V. CONCLUSIONS 



In this paper we have carried out a deeper study of the cosmological behavior of a class 
of modified gravity actions of the form 



\t [bR' + RlBY. 



+ Sm. (74) 



The resulting equations were considerably simplified by choosing as independent variables 
two scalars: the Ricci scalar R and a hnear combination of B? and the Gauss-Bonnet term 
i?QB, namely h R'^ + Rq^. In terms of these two variables, the Einstein equations with matter 
sources take the form given in Eqs. (33, 34). 

The parameter space divides naturally into three regions, b < b, b < b < 0, and b > 
0, where 6 ~ — | (Eq. 60). Although the equations are very stiff, we can nevertheless 
distinguish two kinds of pathological behaviour, for a large region of the parameter space. 
In the region 6 > the presence of a separatrix was shown to be inconsistent with a 
universe whose deceleration parameter q smoothly evolves from 1 at nucleosynthesis to 
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today's negative value. In the region b < b, it was proved that for < g < 1 the universe 
should have hit a spacetime singularity, at which d^Rd^R diverges. This also proves the 
inconsistency of the naive argument which states that ~ i?o implies GR at high redshifts. 

The region b < b < proves to be free of separatrices and singularities. However, a study 
of the propagators of the scalar and tensor modes in backgrounds which start out close to 
GR at nucleosysnthesis, shows that in this interval, for positive integer n, these models are 
affected by classical instabihties. 

In this way we explored the whole parameter space of the generalized models introduced 
in [27], and we found no possibility for a viable cosmology. In [37], the authors explored 
only the simplest case n = 1, and using constraints on the recent expansion history from 
Supernova data, reduced the viable parameter space. In this paper, thanks to the constraints 
on the scalar and tensor propagators introduced in [61], we have eliminated the parameter 
space, leaving no room for these modifications of gravity to be the key component responsible 
for the acceleration of the universe. 

It should be pointed out that the theories we have studied are afflicted by classical 
instabilities (an incorrect sign in the spatial gradients of perturbations), not ghosts (an 
incorrect sign in the time derivatives of perturbations). This means that not even introducing 
a Lorentz-violating cutoff A < 3 MeV along the lines discussed by Cline et al [63] will help. 
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APPENDIX A: EXCEPTIONS 



In this Appendix, we briefly discuss the implications of dropping the assumptions made 
in deriving Eq. (48). 

In principle, one could choose initial conditions at such that the function F2 vanishes 
at that instant, that is equivalent to solve the following algebraic equation for the initial 
condition for A 

^ - ^^i^rnoe-'^2^n^oe-^^-^) ^ ^^^^^ ^ Q ^ ^^^^ 
1 + 7fe. 
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where g is a, function of only b, n, k*. This algebraic equation has at least one and at most 
three real solutions. Therefore, for any fixed value of 6, n, the measure of the subset of the 
solutions which satisfy these initial conditions is zero in the set of all possible solutions, and 
we will not study them any further. 

For 6 < — |, it is impossible that (S^A)^^ vanishes; but it possible, that for a value of k = 
K, ((?kA)k^ — 0. If this happens the previous reasoning discussed in section IV-B regarding 
the impossibility of crossing the singularity, does not apply, as this time oc N — N^. 
However, in this last and this equation gives in turn a relation between b and 

n. In fact one has 

b^b = - ^ [2 + 2n{n + 2) ± (n + 1) V4n^ + 6n + 3] . (A2) 

Having this constraint on b, it is now easier to study the presence of ghosts, classical insta- 
bilities, and superluminal modes only in terms of one parameter n, as ^ is fixed by today's 
data, as already stated in section IV- C. Then, the same result (existence of classical insta- 
bility, i.e. imaginary speed of propagation) applies for b — b{n) , where, also in this case, we 
chose Ni = -10 and 1 < n < 25). 
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